VIS and water analysis
This knitr document contains the code used to analyze and visualize data from the VIS camera system and water data. Knit-ing this document will generate an HTML document that contains both the embedded R code and the output.
Attach the required R packages
#library(qvalue, warn.conflicts = FALSE)
library(ggplot2, warn.conflicts = FALSE)## Warning: package 'ggplot2' was built under R version 3.1.3library(lubridate, warn.conflicts = FALSE)
library(MASS, warn.conflicts = FALSE)
library(car, warn.conflicts = FALSE)## Warning: package 'car' was built under R version 3.1.2library(nlme, warn.conflicts = FALSE)
library(mvtnorm, warn.conflicts = FALSE)
library(grid, warn.conflicts = FALSE)VIS analysis
Download the complete VIS data set. There were a total of 6,399 snapshots with VIS image data, but the download only includes the 6,207 snapshots that were successfully processed by PlantCV. Failed snapshots generally are those that lack a plant (empty pot controls, dead plants, etc.)
#if (!file.exists('vis_snapshots.txt')) {
# download.file('http://files.figshare.com/1845362/vis_snapshots.txt',
# 'vis_snapshots.txt')
#}Read data and format for analysis
vis.data = read.table(file="vis_snapshots_nocorrect.csv", sep=",", header=TRUE)Planting date
planting_date = as.POSIXct("2013-11-26")Add water treatment column coded in barcodes.
vis.data$treatment <- NA
vis.data$treatment[grep("AA", vis.data$plant_id)] <- 100
vis.data$treatment[grep("AB", vis.data$plant_id)] <- 0
vis.data$treatment[grep("AC", vis.data$plant_id)] <- 16
vis.data$treatment[grep("AD", vis.data$plant_id)] <- 33
vis.data$treatment[grep("AE", vis.data$plant_id)] <- 66Add plant genotype column coded in barcodes.
vis.data$genotype <- NA
vis.data$genotype[grep("p1", vis.data$plant_id)] <- 'A10'
vis.data$genotype[grep("p2", vis.data$plant_id)] <- 'B100'
vis.data$genotype[grep("r1", vis.data$plant_id)] <- 'R20'
vis.data$genotype[grep("r2", vis.data$plant_id)] <- 'R70'
vis.data$genotype[grep("r3", vis.data$plant_id)] <- 'R98'
vis.data$genotype[grep("r4", vis.data$plant_id)] <- 'R102'
vis.data$genotype[grep("r5", vis.data$plant_id)] <- 'R128'
vis.data$genotype[grep("r6", vis.data$plant_id)] <- 'R133'
vis.data$genotype[grep("r7", vis.data$plant_id)] <- 'R161'
vis.data$genotype[grep("r8", vis.data$plant_id)] <- 'R187'Add genotype x treatment group column
vis.data$group = paste(vis.data$genotype,'-',vis.data$treatment,sep='')Add calendar-time data column using the Unix-time data
vis.data$date = as.POSIXct(vis.data$datetime, origin = "1970-01-01")Calculate days after planting from planting data
vis.data$dap = as.numeric(vis.data$date - planting_date)Soil volume water content
Use linear regression to create a simple model for using water volume to predict soil volume water content.
Download data for soil wet/dry weight and volume water content measurements.
if (!file.exists('soil_weigth_vwc.txt')) {
download.file('http://files.figshare.com/1939954/soil_weigth_vwc.txt',
'soil_weigth_vwc.txt')
}Read the soil volume water content data
vwc.data = read.table(file="soil_weigth_vwc.txt", sep="\t", header=TRUE)Calculate the average soil dry weight
mean(vwc.data$weight_dry)## [1] 72.92138Create a linear model for water volume and volume water content. Water volumes >= 260 appear to have saturated the soil water carrying capacity, so remove them from the model.
vwc.lm = lm(vwc_wet ~ water_vol, data=vwc.data[vwc.data$water_vol < 260,])
summary(vwc.lm)##
## Call:
## lm(formula = vwc_wet ~ water_vol, data = vwc.data[vwc.data$water_vol <
## 260, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5426 -1.6196 0.3976 1.4214 4.0771
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.032863 1.389638 -2.182 0.0391 *
## water_vol 0.233197 0.008145 28.630 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.082 on 24 degrees of freedom
## Multiple R-squared: 0.9716, Adjusted R-squared: 0.9704
## F-statistic: 819.7 on 1 and 24 DF, p-value: < 2.2e-16Plot the model
vwc.plot = ggplot(vwc.data[vwc.data$water_vol < 260,], aes(x=water_vol, y=vwc_wet)) +
geom_point(size=2) +
geom_smooth(method='lm', formula=y~x) +
scale_x_continuous("Water volume (ml)") +
scale_y_continuous("Volume water content (%)") +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))
print(vwc.plot)
Predict volume water contents for the water treatment groups
treatment.water = data.frame(water_vol=c(217, 144.5, 72))
treatment.water$vwc = predict(object = vwc.lm, newdata=treatment.water)
print(treatment.water)## water_vol vwc
## 1 217.0 47.57084
## 2 144.5 30.66407
## 3 72.0 13.75731Convert VIS camera zoom units
LemnaTec VIS camera zoom units range from 1 to 6000, which correspond to 1 to 6X zoom.
zoom.lm = lm(zoom.camera ~ zoom, data=data.frame(zoom=c(1,6000),
zoom.camera=c(1,6)))
summary(zoom.lm)##
## Call:
## lm(formula = zoom.camera ~ zoom, data = data.frame(zoom = c(1,
## 6000), zoom.camera = c(1, 6)))
##
## Residuals:
## ALL 2 residuals are 0: no residual degrees of freedom!
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9991665 NA NA NA
## zoom 0.0008335 NA NA NA
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 1 and 0 DF, p-value: NAVIS zoom correction
In this section we define models that are used to convert area and length between camera zoom levels to a common scale.
Download data for a reference object imaged at different zoom levels.
#if (!file.exists('zoom_calibration_data.txt')) {
# download.file('http://files.figshare.com/1845355/zoom_calibration_data.txt',
# 'zoom_calibration_data.txt')
#}Read zoom calibrartion data
z.data = read.table(file="zoom_calibration_data.txt", sep="\t", header=TRUE)Calculate px per cm
z.data$px_cm = z.data$length_px / z.data$length_cmConvert LemnaTec zoom units to camera zoom units
z.data$zoom.camera = predict(object = zoom.lm, newdata=z.data)
vis.data$zoom = vis.data$sv_zoom
vis.data$sv.zoom.camera = predict(object = zoom.lm, newdata=vis.data)
vis.data$zoom = vis.data$tv_zoom
vis.data$tv.zoom.camera = predict(object = zoom.lm, newdata=vis.data)Zoom correction for area
Fit a variety of regression models to relative object area by zoom level. Test whether an exponential or 2nd order polynomial model fits best (lowest AIC)
Non-linear regression (exponential)
area.coef = coef(nls(log(rel_area) ~ log(a * exp(b * zoom.camera)),
z.data, start = c(a = 1, b = 0.01)))
area.coef = data.frame(a=area.coef[1], b=area.coef[2])
area.nls = nls(rel_area ~ a * exp(b * zoom.camera),
data = z.data, start=c(a=area.coef$a, b=area.coef$b))
summary(area.nls)##
## Formula: rel_area ~ a * exp(b * zoom.camera)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 0.34409 0.01634 21.06 <2e-16 ***
## b 0.88510 0.01206 73.40 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2738 on 33 degrees of freedom
##
## Number of iterations to convergence: 3
## Achieved convergence tolerance: 1.513e-08Non-linear regression (polynomial)
area.pol = lm(rel_area ~ zoom.camera + I(zoom.camera^2), z.data)
summary(area.pol)##
## Call:
## lm(formula = rel_area ~ zoom.camera + I(zoom.camera^2), data = z.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.7458 -0.4521 0.1182 0.3087 1.6483
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.85322 0.54151 8.962 3.08e-10 ***
## zoom.camera -5.11648 0.45312 -11.292 1.07e-12 ***
## I(zoom.camera^2) 1.73120 0.08478 20.419 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.496 on 32 degrees of freedom
## Multiple R-squared: 0.9892, Adjusted R-squared: 0.9885
## F-statistic: 1465 on 2 and 32 DF, p-value: < 2.2e-16AIC
AIC(area.nls, area.pol)## df AIC
## area.nls 3 12.58188
## area.pol 4 55.10572The exponential model minimizes AIC. Plot exponential model.
nls.plot = ggplot(z.data, aes(x=zoom.camera, y=rel_area)) +
geom_point(size=2.5) +
scale_x_continuous(lim=c(0.5,4.5), "Camera zoom setting") +
scale_y_continuous(lim=c(0,17),
"Reference object relative pixel area") +
stat_smooth(data=z.data, aes(x=zoom.camera, y=rel_area),
method="nls", se=FALSE, formula=y ~ a * exp(b * x),
start=c(a=area.coef$a, b=area.coef$b)) +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))nls.plot = nls.plot + labs(title='Figure 12A')
print(nls.plot)
Zoom correction for length
Fit a variety of regression models to px/cm by zoom level. Test whether an exponential or 2nd order polynomial model fits best (lowest AIC).
Non-linear regression (exponential)
len.coef = coef(nls(log(px_cm) ~ log(a * exp(b * zoom.camera)),
z.data[z.data$camera == 'VIS SV',], start = c(a = 1, b = 0.01)))
len.coef = data.frame(a=len.coef[1], b=len.coef[2])
len.nls = nls(px_cm ~ a * exp(b * zoom.camera),
data = z.data[z.data$camera == 'VIS SV',],
start=c(a=len.coef$a, b=len.coef$b))
summary(len.nls)##
## Formula: px_cm ~ a * exp(b * zoom.camera)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 16.580319 0.271317 61.11 <2e-16 ***
## b 0.449384 0.004676 96.10 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.035 on 15 degrees of freedom
##
## Number of iterations to convergence: 3
## Achieved convergence tolerance: 1.958e-06Length zoom correction using a 2 order polynomial (px/cm given a zoom setting). Length correction only works for side-view images right now because scale in top-view images are affected by the plant lifter position in addition to zoom.
len.poly = lm(px_cm ~ zoom.camera + I(zoom.camera^2),
data=z.data[z.data$camera == 'VIS SV',])
summary(len.poly)##
## Call:
## lm(formula = px_cm ~ zoom.camera + I(zoom.camera^2), data = z.data[z.data$camera ==
## "VIS SV", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.75431 -0.44016 -0.09013 0.40661 1.28788
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.1697 0.9308 33.49 9.13e-15 ***
## zoom.camera -8.9629 0.7901 -11.34 1.92e-08 ***
## I(zoom.camera^2) 6.5778 0.1503 43.75 2.24e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5886 on 14 degrees of freedom
## Multiple R-squared: 0.9996, Adjusted R-squared: 0.9995
## F-statistic: 1.718e+04 on 2 and 14 DF, p-value: < 2.2e-16AIC
AIC(len.nls, len.poly)## df AIC
## len.nls 3 53.28426
## len.poly 4 34.92045The polynomial model minimizes AIC. Plot polynomial model.
poly.plot = ggplot(z.data[z.data$camera == 'VIS SV',], aes(x=zoom.camera, y=px_cm)) +
geom_point(size=2.5) +
scale_x_continuous(lim=c(0.5,4.5), "Camera zoom setting") +
scale_y_continuous(lim=c(0,150),
"Reference object length scale (px/cm)") +
stat_smooth(data=z.data[z.data$camera == 'VIS SV',],
aes(x=zoom.camera, y=px_cm), method="lm",
formula=y ~ x + I(x^2)) +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))poly.plot = poly.plot + labs(title='Figure 12B')
print(poly.plot)
Correct VIS data for differences in camera zoom
Create zoom-corrected VIS data frame
vis.data.zoom = vis.data[,c('plant_id', 'datetime', 'treatment', 'genotype', 'group', 'date', 'dap', 'solidity', 'outlier')]Predict relative area zoom correction factors
vis.data$zoom.camera = vis.data$sv.zoom.camera
vis.data$sv_rel_area = predict(object = area.nls, newdata = vis.data)
vis.data$zoom.camera = vis.data$tv.zoom.camera
vis.data$tv_rel_area = predict(object = area.nls, newdata = vis.data)Calculate total zoom-corrected side-view and top-view area
vis.data.zoom$sv_area = (vis.data$sv0_area / vis.data$sv_rel_area) + (vis.data$sv90_area / vis.data$sv_rel_area) +
(vis.data$sv180_area / vis.data$sv_rel_area) + (vis.data$sv270_area / vis.data$sv_rel_area)
vis.data.zoom$tv_area = vis.data$tv_area / vis.data$tv_rel_areaCalculate zoom-corrected lengths
vis.data$zoom.camera = vis.data$sv.zoom.camera
vis.data$px_cm = predict(object = len.poly, newdata=vis.data)
vis.data.zoom$extent_x = vis.data$extent_x / vis.data$px_cm
vis.data.zoom$extent_y = vis.data$extent_y / vis.data$px_cm
vis.data.zoom$height_above_bound = vis.data$height_above_bound / vis.data$px_cmHeight modeling
In this section we measure how well PlantCV estimates plant height.
Download data manually measured plant height data set (n = 173).
#if (!file.exists('height_model_data.txt')) {
# download.file('http://files.figshare.com/1845361/height_model_data.txt',
# 'height_model_data.txt')
#}Read height data.
manual.ht.data = read.table(file="manual_height_samples.csv",sep=",",header=TRUE)Convert human-readable date to Unix time
manual.ht.data$datetime = as.numeric(ymd_hms(manual.ht.data$timestamp, tz="America/Chicago"))Get height data from the VIS data for each manual height sample
ht.data = merge(manual.ht.data, vis.data.zoom, by = c('plant_id', 'datetime'))Convert LemnaTec zoom units to camera zoom units
ht.data$zoom.camera = predict(object = zoom.lm, newdata=ht.data)Height linear model
height.model = lm(manual_height_cm~height_above_bound, ht.data)
summary(height.model)##
## Call:
## lm(formula = manual_height_cm ~ height_above_bound, data = ht.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.99465 -0.23346 0.01155 0.24932 2.21212
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.372426 0.080998 -4.598 7.7e-06 ***
## height_above_bound 1.024824 0.002947 347.694 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6406 on 193 degrees of freedom
## Multiple R-squared: 0.9984, Adjusted R-squared: 0.9984
## F-statistic: 1.209e+05 on 1 and 193 DF, p-value: < 2.2e-16Plot the height linear model.
hlm.plot = ggplot(ht.data, aes(x=height_above_bound, y=manual_height_cm)) +
geom_point(aes(color=factor(round(zoom.camera,1))),size=2.5) +
geom_smooth(method="lm", formula=y~x, color='black') +
scale_x_continuous(lim=c(0,60), "Estimated height (cm)") +
scale_y_continuous(lim=c(0,60), "Manually measured height (cm)") +
theme_bw() +
theme(legend.position=c(0.2,0.75),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Zoom")hlm.plot = hlm.plot + labs(title='Figure 3A')
print(hlm.plot)
Biomass modeling
In this section we model fresh- and dry-weight above ground biomass using image measurements.
Download data manually measured plant biomass data set (n = 41).
#if (!file.exists('biomass_model_data.txt')) {
# download.file('http://files.figshare.com/1845360/biomass_model_data.txt',
# 'biomass_model_data.txt')
#}Read biomass data.
manual.st.data = read.table(file='manual_biomass_samples.csv', sep=",", header=TRUE, stringsAsFactors=FALSE)Get data from the VIS data for each manual biomass sample
st.data = merge(manual.st.data, vis.data.zoom, by = c('plant_id', 'datetime'))Create out-of-frame indicator variable.
st.data$outind = NA
st.data[st.data$outlier == 'True',]$outind = 1
st.data[st.data$outlier == 'False',]$outind = 0Genotype indicator variable.
st.data$group = NA
st.data[st.data$genotype == 'A10',]$group = 0
st.data[st.data$genotype == 'B100',]$group = 1Fresh-weight biomass
Full model. Includes side-view area, top-view area and height.
fw.full = lm(fresh_weight ~ sv_area * tv_area * height_above_bound, st.data)Step-wise model selection with AIC.
fw.step = stepAIC(fw.full, direction="both")## Start: AIC=-54.68
## fresh_weight ~ sv_area * tv_area * height_above_bound
##
## Df Sum of Sq RSS AIC
## - sv_area:tv_area:height_above_bound 1 0.031061 7.3442 -56.506
## <none> 7.3132 -54.680
##
## Step: AIC=-56.51
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:tv_area +
## sv_area:height_above_bound + tv_area:height_above_bound
##
## Df Sum of Sq RSS AIC
## - sv_area:tv_area 1 0.007503 7.3517 -58.464
## - tv_area:height_above_bound 1 0.059093 7.4033 -58.177
## - sv_area:height_above_bound 1 0.143762 7.4880 -57.711
## <none> 7.3442 -56.506
## + sv_area:tv_area:height_above_bound 1 0.031061 7.3132 -54.680
##
## Step: AIC=-58.46
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:height_above_bound +
## tv_area:height_above_bound
##
## Df Sum of Sq RSS AIC
## - tv_area:height_above_bound 1 0.23599 7.5877 -59.169
## - sv_area:height_above_bound 1 0.31978 7.6715 -58.718
## <none> 7.3517 -58.464
## + sv_area:tv_area 1 0.00750 7.3442 -56.506
##
## Step: AIC=-59.17
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:height_above_bound
##
## Df Sum of Sq RSS AIC
## - sv_area:height_above_bound 1 0.30538 7.8931 -59.551
## <none> 7.5877 -59.169
## + tv_area:height_above_bound 1 0.23599 7.3517 -58.464
## - tv_area 1 0.56816 8.1559 -58.208
## + sv_area:tv_area 1 0.18440 7.4033 -58.177
##
## Step: AIC=-59.55
## fresh_weight ~ sv_area + tv_area + height_above_bound
##
## Df Sum of Sq RSS AIC
## - tv_area 1 0.2631 8.156 -60.207
## <none> 7.893 -59.551
## + sv_area:height_above_bound 1 0.3054 7.588 -59.169
## + tv_area:height_above_bound 1 0.2216 7.671 -58.718
## + sv_area:tv_area 1 0.0478 7.845 -57.800
## - height_above_bound 1 0.8110 8.704 -57.541
## - sv_area 1 23.8156 31.709 -4.536
##
## Step: AIC=-60.21
## fresh_weight ~ sv_area + height_above_bound
##
## Df Sum of Sq RSS AIC
## <none> 8.156 -60.207
## + tv_area 1 0.263 7.893 -59.551
## - height_above_bound 1 0.654 8.810 -59.046
## + sv_area:height_above_bound 1 0.000 8.156 -58.208
## - sv_area 1 157.554 165.710 61.263summary(fw.step)##
## Call:
## lm(formula = fresh_weight ~ sv_area + height_above_bound, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.95619 -0.19476 -0.00613 0.06206 1.49454
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.455e-02 2.002e-01 0.173 0.864
## sv_area 3.972e-05 1.466e-06 27.093 <2e-16 ***
## height_above_bound -4.029e-02 2.309e-02 -1.745 0.089 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4633 on 38 degrees of freedom
## Multiple R-squared: 0.9842, Adjusted R-squared: 0.9833
## F-statistic: 1181 on 2 and 38 DF, p-value: < 2.2e-16AIC model
fw.aic = lm(fresh_weight ~ sv_area + tv_area + height_above_bound +
sv_area*height_above_bound, st.data)
summary(fw.aic)##
## Call:
## lm(formula = fresh_weight ~ sv_area + tv_area + height_above_bound +
## sv_area * height_above_bound, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.14016 -0.10203 -0.00920 0.07919 1.56414
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.812e-01 2.356e-01 0.769 0.4468
## sv_area 4.538e-05 4.288e-06 10.584 1.34e-12 ***
## tv_area -1.026e-05 6.248e-06 -1.642 0.1093
## height_above_bound -6.232e-02 2.709e-02 -2.301 0.0273 *
## sv_area:height_above_bound 1.924e-07 1.598e-07 1.204 0.2366
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4591 on 36 degrees of freedom
## Multiple R-squared: 0.9853, Adjusted R-squared: 0.9836
## F-statistic: 602 on 4 and 36 DF, p-value: < 2.2e-16The AIC model contains tv_area and height which does not have a significant coefficient, test dropping.
fw.red = lm(fresh_weight ~ sv_area, st.data)
summary(fw.red)##
## Call:
## lm(formula = fresh_weight ~ sv_area, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.10964 -0.16717 0.00494 0.24429 1.30677
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.704e-01 1.002e-01 -2.70 0.0102 *
## sv_area 3.755e-05 7.931e-07 47.34 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4753 on 39 degrees of freedom
## Multiple R-squared: 0.9829, Adjusted R-squared: 0.9825
## F-statistic: 2241 on 1 and 39 DF, p-value: < 2.2e-16Goodness of fit.
anova(fw.aic, fw.red)## Analysis of Variance Table
##
## Model 1: fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area *
## height_above_bound
## Model 2: fresh_weight ~ sv_area
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 36 7.5877
## 2 39 8.8099 -3 -1.2222 1.9329 0.1417SV area model.
sv.model = lm(fresh_weight ~ sv_area, st.data)
summary(sv.model)##
## Call:
## lm(formula = fresh_weight ~ sv_area, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.10964 -0.16717 0.00494 0.24429 1.30677
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.704e-01 1.002e-01 -2.70 0.0102 *
## sv_area 3.755e-05 7.931e-07 47.34 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4753 on 39 degrees of freedom
## Multiple R-squared: 0.9829, Adjusted R-squared: 0.9825
## F-statistic: 2241 on 1 and 39 DF, p-value: < 2.2e-16Plot SV model
sv.model.plot = ggplot(st.data,aes(x=sv_area/1e5, y=fresh_weight)) +
geom_smooth(method="lm", color="black", formula = y ~ x) +
geom_point(size=2.5) +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Fresh-weight biomass (g)") +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))sv.model.plot = sv.model.plot + labs(title='Figure 4A')
print(sv.model.plot)
SV area model with out-of-frame.
sv.ind.model = lm(fresh_weight ~ sv_area + outind, st.data)
anova(sv.model, sv.ind.model)## Analysis of Variance Table
##
## Model 1: fresh_weight ~ sv_area
## Model 2: fresh_weight ~ sv_area + outind
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 39 8.8099
## 2 38 7.0902 1 1.7196 9.2164 0.004316 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(sv.ind.model)##
## Call:
## lm(formula = fresh_weight ~ sv_area + outind, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.78369 -0.19699 -0.01344 0.21652 1.20823
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.202e-01 9.252e-02 -2.381 0.02241 *
## sv_area 3.831e-05 7.636e-07 50.174 < 2e-16 ***
## outind -5.242e-01 1.727e-01 -3.036 0.00432 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.432 on 38 degrees of freedom
## Multiple R-squared: 0.9862, Adjusted R-squared: 0.9855
## F-statistic: 1361 on 2 and 38 DF, p-value: < 2.2e-16Plot SV model with out-of-frame.
sv.model.out.plot = ggplot(st.data, aes(x=sv_area/1e5, y=fresh_weight,
group=outlier, color=outlier)) +
geom_point(size=2.5) +
geom_smooth(method="lm", color="black") +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Fresh-weight biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Out-of-frame")
print(sv.model.out.plot)
SV area model with genotype
sv.gt.model = lm(fresh_weight ~ sv_area + group, st.data)
summary(sv.gt.model)##
## Call:
## lm(formula = fresh_weight ~ sv_area + group, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.24908 -0.29687 0.01704 0.21582 1.17723
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.150e-01 1.266e-01 -3.279 0.00223 **
## sv_area 3.775e-05 7.797e-07 48.414 < 2e-16 ***
## group 2.614e-01 1.460e-01 1.790 0.08135 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4624 on 38 degrees of freedom
## Multiple R-squared: 0.9842, Adjusted R-squared: 0.9834
## F-statistic: 1186 on 2 and 38 DF, p-value: < 2.2e-16Plot SV model with genotype.
sv.model.gen.plot = ggplot(st.data, aes(x=sv_area/1e5, y=fresh_weight,
group=genotype, color=genotype)) +
geom_point(size=2.5) +
geom_smooth(method="lm", color="black") +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Fresh-weight biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype")
print(sv.model.gen.plot)
Dry-weight biomass
SV area model
dry.sv.model = lm(dry_weight ~ sv_area, st.data)
summary(dry.sv.model)##
## Call:
## lm(formula = dry_weight ~ sv_area, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16009 -0.04743 0.01569 0.04033 0.18561
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.915e-02 1.389e-02 -3.539 0.00106 **
## sv_area 4.425e-06 1.100e-07 40.233 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06591 on 39 degrees of freedom
## Multiple R-squared: 0.9765, Adjusted R-squared: 0.9759
## F-statistic: 1619 on 1 and 39 DF, p-value: < 2.2e-16Plot dry-weight side-view model
dry.sv.model.plot = ggplot(st.data,aes(x=sv_area/1e5, y=dry_weight)) +
geom_smooth(method="lm", color="black", formula = y ~ x) +
geom_point(size=2.5) +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Dry-weight biomass (g)") +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))dry.sv.model.plot = dry.sv.model.plot + labs(title='Supplemental Figure S2')
print(dry.sv.model.plot)
SV area model with out-of-frame.
dry.sv.ind.model = lm(dry_weight ~ sv_area + outind, st.data)
summary(dry.sv.ind.model)##
## Call:
## lm(formula = dry_weight ~ sv_area + outind, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.112715 -0.030983 0.006637 0.024625 0.127446
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.052e-02 1.196e-02 -3.388 0.001650 **
## sv_area 4.556e-06 9.869e-08 46.166 < 2e-16 ***
## outind -9.023e-02 2.232e-02 -4.043 0.000248 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05583 on 38 degrees of freedom
## Multiple R-squared: 0.9836, Adjusted R-squared: 0.9827
## F-statistic: 1136 on 2 and 38 DF, p-value: < 2.2e-16dry.sv.model.out.plot = ggplot(st.data, aes(x=sv_area/1e5, y=dry_weight,
group=outlier, color=outlier)) +
geom_point(size=2.5) +
geom_smooth(method="lm", color="black") +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Dry-weight biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Out-of-frame")
print(dry.sv.model.out.plot)
SV area model with genotypes
dry.sv.gt.model = lm(dry_weight ~ sv_area + group, st.data)
summary(dry.sv.gt.model)##
## Call:
## lm(formula = dry_weight ~ sv_area + group, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.15945 -0.04653 0.01473 0.03951 0.18609
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.015e-02 1.827e-02 -2.745 0.0092 **
## sv_area 4.426e-06 1.126e-07 39.315 <2e-16 ***
## group 1.813e-03 2.108e-02 0.086 0.9319
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06676 on 38 degrees of freedom
## Multiple R-squared: 0.9765, Adjusted R-squared: 0.9752
## F-statistic: 788.8 on 2 and 38 DF, p-value: < 2.2e-16dry.sv.model.gen.plot = ggplot(st.data, aes(x=sv_area/1e5, y=dry_weight,
group=genotype, color=genotype)) +
geom_point(size=2.5) +
geom_smooth(method="lm", color="black") +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Dry-weight biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype")
print(dry.sv.model.gen.plot)
Analyze VIS data
In this section we analyze geometric traits from the complete VIS data set.
Remove plants that were sampled for biomass measurements.
vis.data.zoom = vis.data.zoom[!vis.data.zoom$plant_id %in% st.data$plant_id,]Predict fresh- and dry-weight biomass from linear models
vis.data.zoom$fw_biomass = predict.lm(object = sv.model, newdata=vis.data.zoom)
vis.data.zoom$dw_biomass = predict.lm(object = dry.sv.model, newdata=vis.data.zoom)Plant height analysis
Plot height for S. viridis and S. italica water treatments 100% and 33% full-capacity.
height.plot = ggplot(vis.data.zoom[(vis.data.zoom$genotype == 'A10' |
vis.data.zoom$genotype == 'B100') &
(vis.data.zoom$treatment == 100 |
vis.data.zoom$treatment == 33),],
aes(x=dap, y=height_above_bound, color=factor(group))) +
geom_point(size=2.5) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(0,60),
name="Estimated height (cm)") +
theme_bw() +
theme(legend.position=c(0.25,0.75),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype-treatment")height.plot = height.plot + labs(title='Figure 3B')
print(height.plot)
Plot height for all genotypes.
height.facet.plot = ggplot(vis.data.zoom[vis.data.zoom$treatment == 100 |
vis.data.zoom$treatment == 33,],
aes(x=dap, y=height_above_bound, color=factor(treatment))) +
geom_point(size=1) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(0,60),
name="Estimated height (cm)") +
theme_bw() +
theme(legend.position=c(0.8,0.1),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Treatment") +
facet_wrap(~ genotype, ncol = 3)height.facet.plot = height.facet.plot + labs(title='Supplemental Figure S1')
print(height.facet.plot)
Genotype differences in height
For each day, calculate the mean and 95% confidence intervals for each genotype in the control water group.
height_per_day = function(vis.data) {
dap = c()
genotype = c()
int.low = c()
int.up = c()
est = c()
# Loop through each day and genotype
for(d in min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))) {
if(d != 24) {
for(g in levels(factor(vis.data$genotype))) {
h.data = vis.data[vis.data$genotype == g & vis.data$treatment == 100 &
as.integer(vis.data$dap) == d,]$height_above_bound
if(sd(h.data) != 0) {
dap = c(dap,d)
genotype = c(genotype,g)
test = t.test(h.data)
int.low = c(int.low,test$conf.int[1])
int.up = c(int.up,test$conf.int[2])
est = c(est,test$estimate)
}
}
}
}
#drought_resp = drought_resp[-14,]
results = data.frame(dap=as.numeric(dap),
genotype=genotype,
conf.int.low=int.low,
conf.int.up=int.up,
mean=est)
return(results)
}Height per day 95% CI. Write results to file.
height.perday.results = height_per_day(vis.data.zoom)
write.csv(height.perday.results,file="height.perday.100water.csv")
print(height.perday.results)## dap genotype conf.int.low conf.int.up mean
## 1 11 A10 4.1138046 5.278839 4.696322
## 2 11 B100 2.5259319 3.326661 2.926297
## 3 11 R102 1.2684481 3.282102 2.275275
## 4 11 R128 1.6661243 2.966349 2.316237
## 5 11 R133 2.8262168 5.442689 4.134453
## 6 11 R161 2.8964352 3.279606 3.088021
## 7 11 R187 2.2374378 4.006971 3.122204
## 8 11 R20 4.2648901 6.073519 5.169204
## 9 11 R70 3.3847556 4.309511 3.847133
## 10 11 R98 5.5427726 6.565899 6.054336
## 11 12 A10 4.6895754 6.770786 5.730181
## 12 12 B100 2.7955379 4.609936 3.702737
## 13 12 R102 2.0129905 5.055558 3.534274
## 14 12 R128 2.3173878 4.551756 3.434572
## 15 12 R133 0.1833279 16.040642 8.111985
## 16 12 R161 2.7075759 5.056310 3.881943
## 17 12 R187 4.3345350 6.973768 5.654151
## 18 12 R20 -12.0782382 28.023729 7.972745
## 19 12 R70 3.0529634 5.786161 4.419562
## 20 12 R98 7.0338485 9.288841 8.161345
## 21 13 A10 7.2633774 8.933993 8.098685
## 22 13 B100 5.4971969 6.483222 5.990209
## 23 13 R102 3.9400526 6.147069 5.043561
## 24 13 R128 4.8671575 6.363675 5.615416
## 25 13 R133 7.5927485 9.939040 8.765894
## 26 13 R161 5.6093371 7.035664 6.322501
## 27 13 R187 5.3732471 7.338795 6.356021
## 28 13 R20 7.2176543 10.572329 8.894992
## 29 13 R70 6.2027656 7.171094 6.686930
## 30 13 R98 8.7080117 11.816909 10.262460
## 31 14 A10 7.2992964 10.650554 8.974925
## 32 14 B100 7.2095969 8.594076 7.901836
## 33 14 R102 4.6619405 8.705043 6.683492
## 34 14 R128 5.6149445 9.147871 7.381408
## 35 14 R133 11.8870482 12.800615 12.343831
## 36 14 R161 5.5132820 8.133326 6.823304
## 37 14 R187 7.7707801 10.776425 9.273602
## 38 14 R20 4.4397218 17.446651 10.943186
## 39 14 R70 6.2590147 10.235183 8.247099
## 40 14 R98 11.0626447 12.903724 11.983184
## 41 15 A10 10.5250697 12.024864 11.274967
## 42 15 B100 8.6034387 10.159156 9.381297
## 43 15 R102 5.5131064 8.935807 7.224457
## 44 15 R128 8.3515328 10.617965 9.484749
## 45 15 R133 13.4816282 16.926665 15.204147
## 46 15 R161 8.5636321 10.277182 9.420407
## 47 15 R187 7.8971054 11.444819 9.670962
## 48 15 R20 12.5540041 16.579996 14.567000
## 49 15 R70 9.4534249 11.399830 10.426628
## 50 15 R98 14.4840908 15.932854 15.208472
## 51 16 A10 10.3183950 13.664952 11.991674
## 52 16 B100 10.3217123 12.069260 11.195486
## 53 16 R102 7.4126348 10.542533 8.977584
## 54 16 R128 8.2558336 15.124237 11.690035
## 55 16 R133 16.4716388 17.833934 17.152786
## 56 16 R161 8.2040259 11.746160 9.975093
## 57 16 R187 11.1664495 13.991223 12.578836
## 58 16 R20 4.4901376 25.384677 14.937407
## 59 16 R70 10.1166241 15.935022 13.025823
## 60 16 R98 16.3400141 17.763702 17.051858
## 61 17 A10 13.8288699 16.273493 15.051181
## 62 17 B100 12.8021984 14.751262 13.776730
## 63 17 R102 4.6784026 20.775747 12.727075
## 64 17 R128 13.3720766 17.669471 15.520774
## 65 17 R133 17.4574601 17.473869 17.465665
## 66 17 R161 9.2299092 20.579302 14.904606
## 67 17 R187 13.3249948 15.732300 14.528647
## 68 17 R20 17.2492825 17.540747 17.395015
## 69 17 R70 13.8775205 16.001331 14.939426
## 70 17 R98 17.4389051 17.482331 17.460618
## 71 18 A10 14.4288846 16.848835 15.638860
## 72 18 B100 15.1035778 16.693924 15.898751
## 73 18 R102 11.9392515 17.839681 14.889466
## 74 18 R128 14.4065180 16.592639 15.499579
## 75 18 R161 11.3654767 16.590256 13.977866
## 76 18 R187 14.2918010 18.966135 16.628968
## 77 18 R20 17.3669600 17.559323 17.463141
## 78 18 R70 16.2359873 18.115003 17.175495
## 79 18 R98 17.4695785 17.486983 17.478281
## 80 19 A10 15.5274345 17.339111 16.433273
## 81 19 B100 14.1045628 16.056801 15.080682
## 82 19 R102 9.2740593 18.382349 13.828204
## 83 19 R128 17.4682044 17.481868 17.475036
## 84 19 R133 17.4691725 17.478305 17.473739
## 85 19 R161 13.7376679 17.711630 15.724649
## 86 19 R187 15.1000979 17.975410 16.537754
## 87 19 R20 17.4722809 17.482904 17.477592
## 88 19 R70 16.3674982 17.915366 17.141432
## 89 20 A10 15.4815153 17.904691 16.693103
## 90 20 B100 14.9940006 18.062583 16.528292
## 91 20 R102 14.8917746 18.071450 16.481612
## 92 20 R128 15.8643804 18.103082 16.983731
## 93 20 R133 17.4653321 17.482145 17.473739
## 94 20 R161 13.7669020 17.208033 15.487467
## 95 20 R187 16.8032817 17.759659 17.281470
## 96 20 R20 17.4540444 17.497470 17.475757
## 97 20 R70 15.5391179 18.384945 16.962032
## 98 20 R98 17.4675530 17.483962 17.475757
## 99 21 A10 18.0979070 22.866964 20.482435
## 100 21 B100 17.2659492 20.740679 19.003314
## 101 21 R102 10.3649359 27.239881 18.802408
## 102 21 R128 23.2259512 29.985433 26.605692
## 103 21 R133 35.2919365 38.324288 36.808112
## 104 21 R161 15.6427963 23.752576 19.697686
## 105 21 R187 20.9302871 28.030502 24.480395
## 106 21 R20 28.5799826 34.970515 31.775249
## 107 21 R70 22.9352632 26.500120 24.717692
## 108 21 R98 28.6623854 35.443911 32.053148
## 109 22 A10 20.8063040 25.919929 23.363117
## 110 22 B100 17.8425711 25.156676 21.499623
## 111 22 R102 16.7256650 28.193792 22.459728
## 112 22 R128 22.0557320 29.847211 25.951471
## 113 22 R133 33.3894751 49.007692 41.198583
## 114 22 R161 16.4579791 20.578905 18.518442
## 115 22 R187 25.3813857 35.138840 30.260113
## 116 22 R20 23.4710805 58.736805 41.103943
## 117 22 R70 19.0935859 36.220406 27.656996
## 118 22 R98 33.3321134 40.188754 36.760433
## 119 23 A10 23.7887258 29.592135 26.690430
## 120 23 B100 20.1451618 25.383980 22.764571
## 121 23 R102 11.4174069 32.037610 21.727508
## 122 23 R128 29.2037138 35.281459 32.242586
## 123 23 R133 41.3614576 45.130164 43.245811
## 124 23 R161 20.4242532 25.227423 22.825838
## 125 23 R187 23.1617592 38.261040 30.711400
## 126 23 R20 34.0321909 44.150975 39.091583
## 127 23 R70 30.2497912 34.894308 32.572050
## 128 23 R98 31.3077114 41.034772 36.171242
## 129 25 A10 28.8418538 34.316538 31.579196
## 130 25 B100 21.7374532 31.214009 26.475731
## 131 25 R102 19.3852591 38.594612 28.989935
## 132 25 R128 30.4325931 40.029886 35.231239
## 133 25 R133 47.2030660 54.937143 51.070105
## 134 25 R161 21.5395173 26.577798 24.058658
## 135 25 R187 35.8052978 41.428971 38.617134
## 136 25 R20 24.4029064 70.696036 47.549471
## 137 25 R70 28.2682884 48.061912 38.165100
## 138 25 R98 40.7567043 49.171868 44.964286
## 139 26 A10 31.3246339 38.418576 34.871605
## 140 26 B100 26.7993073 34.574345 30.686826
## 141 26 R102 14.9726330 46.784895 30.878764
## 142 26 R128 35.9731220 44.331988 40.152555
## 143 26 R133 52.6693551 52.719093 52.694224
## 144 26 R161 23.0900303 39.985494 31.537762
## 145 26 R187 33.1394226 43.142960 38.141191
## 146 26 R20 49.2489502 53.411285 51.330117
## 147 26 R70 41.7251538 46.738979 44.232066
## 148 26 R98 34.7437085 49.954721 42.349215
## 149 27 A10 35.1602826 41.245993 38.203138
## 150 27 B100 27.0240390 37.926865 32.475452
## 151 27 R102 23.9486786 47.091606 35.520142
## 152 27 R128 35.9592946 47.623371 41.791333
## 153 27 R133 52.6385536 52.932629 52.785592
## 154 27 R161 24.2442800 30.727013 27.485647
## 155 27 R187 41.2496516 44.943105 43.096378
## 156 27 R20 48.8256683 55.318931 52.072300
## 157 27 R70 34.9404743 52.806382 43.873428
## 158 27 R98 47.5987571 53.069042 50.333900
## 159 28 A10 37.3444960 43.747626 40.546061
## 160 28 B100 31.0852941 39.304128 35.194711
## 161 28 R102 19.7785402 51.925226 35.851883
## 162 28 R128 45.7837167 50.042465 47.913091
## 163 28 R133 52.6330116 52.823749 52.728380
## 164 28 R161 28.0109976 44.448898 36.229948
## 165 28 R187 35.2088327 47.003038 41.105935
## 166 28 R20 52.5767078 52.970566 52.773637
## 167 28 R70 49.9879362 52.206923 51.097429
## 168 28 R98 42.0724187 54.742865 48.407642
## 169 29 A10 40.2283979 45.252515 42.740457
## 170 29 B100 33.4612448 46.507641 39.984443
## 171 29 R102 26.5575297 53.548336 40.052933
## 172 29 R128 40.3678421 53.396169 46.882005
## 173 29 R133 52.6874775 53.051070 52.869274
## 174 29 R161 31.3730339 37.445685 34.409360
## 175 29 R187 44.0292757 47.520641 45.774958
## 176 29 R20 52.6909028 52.987872 52.839387
## 177 29 R70 42.0424398 56.607023 49.324731
## 178 29 R98 52.0544743 52.857213 52.455843
## 179 30 A10 41.1590692 48.054219 44.606644
## 180 30 B100 36.0961707 44.787331 40.441751
## 181 30 R102 27.1510111 54.451174 40.801093
## 182 30 R128 51.7339576 52.709507 52.221732
## 183 30 R133 52.6812412 52.877987 52.779614
## 184 30 R161 37.2649692 46.088567 41.676768
## 185 30 R187 39.0026350 49.389770 44.196202
## 186 30 R20 52.7731212 52.959449 52.866285
## 187 30 R70 52.6255605 52.942207 52.783884
## 188 30 R98 48.3607663 54.312277 51.336522
## 189 31 A10 44.9409479 48.995137 46.968042
## 190 31 B100 39.6784848 50.930006 45.304246
## 191 31 R102 32.3953484 54.156046 43.275697
## 192 31 R128 45.9896516 54.628337 50.308994
## 193 31 R133 52.8090833 53.001192 52.905138
## 194 31 R161 35.3843979 41.461346 38.422872
## 195 31 R187 46.3300941 49.725217 48.027656
## 196 31 R20 52.5361938 53.321900 52.929047
## 197 31 R70 45.1991354 56.510708 50.854922
## 198 31 R98 52.7570164 52.941683 52.849349
## 199 32 A10 43.9442034 51.699444 47.821823
## 200 32 B100 42.0255866 49.295184 45.660385
## 201 32 R102 33.9379036 55.911966 44.924935
## 202 32 R128 53.2778901 54.523848 53.900869
## 203 32 R133 54.3529911 54.469816 54.411403
## 204 32 R161 39.2900732 49.085036 44.187554
## 205 32 R187 43.5507291 49.570772 46.560751
## 206 32 R20 54.4236761 54.468606 54.446141
## 207 32 R70 54.3062174 54.417340 54.361779
## 208 32 R98 53.3227703 54.787243 54.055007
## 209 33 A10 48.7441481 52.620165 50.682157
## 210 33 B100 45.2839248 53.441351 49.362638
## 211 33 R102 38.4873097 55.935265 47.211287
## 212 33 R128 48.7761391 55.972918 52.374528
## 213 33 R133 54.3633798 54.459427 54.411403
## 214 33 R161 41.5367102 45.138371 43.337541
## 215 33 R187 47.6165089 50.414019 49.015264
## 216 33 R20 54.3071074 54.578858 54.442983
## 217 33 R70 46.8503717 57.968160 52.409266
## 218 33 R98 54.3418212 54.470459 54.406140Treatment differences in height
Statistical analysis for height differences. Use pairwise comparisons on S. viridis and S. italica over two-day intervals.
analyze_height = function(vis.data, genotype) {
days = c()
diff.low = c()
diff.up = c()
pvals = c()
for(day in levels(factor(as.integer(vis.data$dap)))) {
day = as.integer(day)
control = vis.data[(as.integer(vis.data$dap) == day |
as.integer(vis.data$dap) == day + 1) &
vis.data$genotype == genotype &
vis.data$treatment == 100,]$height_above_bound
drought = vis.data[(as.integer(vis.data$dap) == day |
as.integer(vis.data$dap) == day + 1) &
vis.data$genotype == genotype &
vis.data$treatment == 33,]$height_above_bound
test = t.test(x=control, y=drought)
days = c(days, day)
diff.low = c(diff.low, test$conf.int[1])
diff.up = c(diff.up, test$conf.int[2])
pvals = c(pvals, test$p.value)
}
results = data.frame(dap=as.numeric(days),
conf.int.low=diff.low,
conf.int.up=diff.up,
pvalue=pvals)
return(results)
}Test for treatment effect on two-day intervals
a10.height.results = analyze_height(vis.data.zoom, 'A10')
b100.height.results = analyze_height(vis.data.zoom, 'B100')Control for multiple testing by controlling the FDR
qvalues.height = p.adjust(c(a10.height.results$pvalue, b100.height.results$pvalue),method="fdr")
a10.height.results$qvalue = qvalues.height[1:nrow(a10.height.results)]
b100.height.results$qvalue = qvalues.height[
(nrow(a10.height.results)+1):(nrow(a10.height.results) + nrow(b100.height.results))]Assign genotypes for merged table
a10.height.results$genotype = 'A10'
b100.height.results$genotype = 'B100'
print(a10.height.results)## dap conf.int.low conf.int.up pvalue qvalue genotype
## 1 11 0.77273960 2.643169 9.582125e-04 0.0032431807 A10
## 2 12 -0.29914242 2.422404 1.206740e-01 0.2097086778 A10
## 3 13 -0.08842057 2.115835 7.039538e-02 0.1534676662 A10
## 4 14 -0.32917060 2.580711 1.239188e-01 0.2097086778 A10
## 5 15 0.28426142 2.919397 1.924394e-02 0.0529208463 A10
## 6 16 0.14666528 4.378048 3.748556e-02 0.0868086644 A10
## 7 17 -0.87910116 3.018279 2.505424e-01 0.3674622138 A10
## 8 18 0.48991495 3.897893 1.460564e-02 0.0428432116 A10
## 9 19 1.15428894 4.393747 1.873909e-03 0.0058894273 A10
## 10 20 2.29234606 6.796709 2.596630e-04 0.0009520977 A10
## 11 21 3.42568864 8.778005 6.793406e-05 0.0003736373 A10
## 12 22 5.05021107 10.781276 3.072841e-06 0.0001352050 A10
## 13 23 4.77919478 12.771985 2.155551e-04 0.0008622202 A10
## 14 25 3.75450813 10.565738 1.699507e-04 0.0008308701 A10
## 15 26 4.41239580 12.298773 1.911306e-04 0.0008409746 A10
## 16 27 5.94525469 13.501681 1.946187e-05 0.0001951542 A10
## 17 28 5.54257526 13.398589 5.606520e-05 0.0003736373 A10
## 18 29 5.90446577 12.852781 7.632409e-06 0.0001679130 A10
## 19 30 6.17854616 13.985667 2.217662e-05 0.0001951542 A10
## 20 31 6.60003839 14.577970 1.300773e-05 0.0001907800 A10
## 21 32 6.21299029 15.161443 6.320013e-05 0.0003736373 A10
## 22 33 1.73822539 19.572166 2.690657e-02 0.0696405324 A10print(b100.height.results)## dap conf.int.low conf.int.up pvalue qvalue genotype
## 1 11 -1.2594720 0.06992054 0.07673383 0.15346767 B100
## 2 12 -0.9886257 1.13640045 0.88658293 0.88658293 B100
## 3 13 -1.8408812 0.31467086 0.15621418 0.24547943 B100
## 4 14 -1.1123323 1.64762375 0.68794161 0.70394026 B100
## 5 15 -2.6455439 0.31836058 0.11604493 0.20970868 B100
## 6 16 -0.7864901 2.70009313 0.26085877 0.37025116 B100
## 7 17 -2.5824550 1.44347769 0.54836630 0.60320293 B100
## 8 18 -1.0746677 2.55761958 0.40054959 0.51835830 B100
## 9 19 -1.3873902 2.36907573 0.58999371 0.63316398 B100
## 10 20 -0.1953884 3.88187425 0.07458440 0.15346767 B100
## 11 21 -4.7055442 2.44925827 0.51436328 0.60320293 B100
## 12 22 -4.1699221 2.69196199 0.66036184 0.69180764 B100
## 13 23 -5.6805044 2.80840878 0.46303807 0.56593541 B100
## 14 25 -6.9139710 1.33553502 0.17774764 0.26968607 B100
## 15 26 -8.9304594 -0.52709508 0.02867192 0.07008691 B100
## 16 27 -7.3841917 1.12487918 0.14381616 0.23436707 B100
## 17 28 -8.7690652 0.63533080 0.08754828 0.16748367 B100
## 18 29 -6.1716464 2.70162092 0.43167803 0.54268095 B100
## 19 30 -6.3500065 2.33465144 0.35381611 0.48649716 B100
## 20 31 -2.6645326 4.95701544 0.54451691 0.60320293 B100
## 21 32 -2.8011143 5.16714436 0.54834969 0.60320293 B100
## 22 33 -3.1147374 7.53586129 0.38830762 0.51774350 B100Output the A10 and B100 tables to the same file
write.table(a10.height.results, file='height.stats.csv', sep = ',',
row.names = FALSE, append = FALSE)
write.table(b100.height.results, file='height.stats.csv', sep = ',',
row.names = FALSE, append = TRUE, col.names = FALSE)Biomass analysis
Plot fresh-weight biomass for S. viridis and S. italica water treatments 100% and 33% full-capacity.
biomass.plot = ggplot(vis.data.zoom[(vis.data.zoom$genotype == 'A10' |
vis.data.zoom$genotype == 'B100') &
(vis.data.zoom$treatment == 100 |
vis.data.zoom$treatment == 33),],
aes(x=dap, y=fw_biomass, color=factor(group))) +
geom_point(size=2.5) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-1,41),
name="Estimated biomass (g)") +
theme_bw() +
theme(legend.position=c(0.25,0.75),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype-treatment")biomass.plot = biomass.plot + labs(title='Figure 4B')
print(biomass.plot)
Plot fresh-weight biomass for all genotypes.
biomass.facet.plot = ggplot(vis.data.zoom[vis.data.zoom$treatment == 100 |
vis.data.zoom$treatment == 33,],
aes(x=dap, y=fw_biomass, color=factor(treatment))) +
geom_point(size=1) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-1,41),
name="Estimated biomass (g)") +
theme_bw() +
theme(legend.position=c(0.8,0.1),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Treatment") +
facet_wrap(~ genotype, ncol = 3)biomass.facet.plot = biomass.facet.plot + labs(title='Supplemental Figure S3')
print(biomass.facet.plot)
Treatment differences in biomass
Statistical analysis for fresh-weight biomass differences. Use pairwise comparisons on S. viridis and S. italica over two-day intervals.
analyze_biomass = function(vis.data, genotype) {
days = c()
diff.low = c()
diff.up = c()
pvals = c()
for(day in levels(factor(as.integer(vis.data$dap)))) {
day = as.integer(day)
control = vis.data[(as.integer(vis.data$dap) == day |
as.integer(vis.data$dap) == day + 1) &
vis.data$genotype == genotype &
vis.data$treatment == 100,]$fw_biomass
drought = vis.data[(as.integer(vis.data$dap) == day |
as.integer(vis.data$dap) == day + 1) &
vis.data$genotype == genotype &
vis.data$treatment == 33,]$fw_biomass
# If the control group has a lot more replicates (e.g. A10 and B100), randomly sample to drought sample size
if (genotype == 'A10' | genotype == 'B100') {
control = sample(control, size = length(drought))
}
test = t.test(x=control, y=drought)
days = c(days, day)
diff.low = c(diff.low, test$conf.int[1])
diff.up = c(diff.up, test$conf.int[2])
pvals = c(pvals, test$p.value)
}
results = data.frame(dap=as.numeric(days),
conf.int.low=diff.low,
conf.int.up=diff.up,
pvalue=pvals)
return(results)
}Test for treatment effect on two-day intervals
a10.biomass.results = analyze_biomass(vis.data.zoom, 'A10')
b100.biomass.results = analyze_biomass(vis.data.zoom, 'B100')Control for multiple testing by controlling the FDR
qvalues.biomass = p.adjust(c(a10.biomass.results$pvalue, b100.biomass.results$pvalue),method="fdr")
a10.biomass.results$qvalue = qvalues.biomass[1:nrow(a10.biomass.results)]
b100.biomass.results$qvalue = qvalues.biomass[
(nrow(a10.biomass.results)+1):(nrow(a10.biomass.results) +
nrow(b100.biomass.results))]Assign genotypes for merged table
a10.biomass.results$genotype = 'A10'
b100.biomass.results$genotype = 'B100'
print(a10.biomass.results)## dap conf.int.low conf.int.up pvalue qvalue genotype
## 1 11 0.06618061 0.2366926 1.407098e-03 2.948206e-03 A10
## 2 12 -0.08602355 0.2214458 3.733164e-01 4.211774e-01 A10
## 3 13 -0.14878889 0.2530395 6.000685e-01 6.286432e-01 A10
## 4 14 -0.11355314 0.4706694 2.189056e-01 2.832896e-01 A10
## 5 15 -0.16532964 0.7004685 2.138495e-01 2.832896e-01 A10
## 6 16 -0.13466603 1.1361927 1.146885e-01 1.802249e-01 A10
## 7 17 -0.71306484 2.5798564 2.404504e-01 3.022805e-01 A10
## 8 18 1.01265567 2.7833441 1.984776e-04 5.822009e-04 A10
## 9 19 0.59986296 2.8235229 3.885584e-03 7.433292e-03 A10
## 10 20 1.51836098 4.7372845 6.640184e-04 1.623156e-03 A10
## 11 21 2.11154234 6.0531431 3.243150e-04 8.918661e-04 A10
## 12 22 4.28169719 8.8377380 9.412210e-06 3.185671e-05 A10
## 13 23 4.94273400 10.4089661 6.069847e-05 1.907666e-04 A10
## 14 25 6.47172203 11.7791467 6.402658e-07 3.130189e-06 A10
## 15 26 7.10409340 11.8459095 3.113334e-08 1.712334e-07 A10
## 16 27 8.85914954 12.3725610 6.213467e-12 4.556542e-11 A10
## 17 28 9.61069065 12.9472564 5.285230e-13 5.813753e-12 A10
## 18 29 8.91289326 12.3369483 3.464972e-12 3.049175e-11 A10
## 19 30 10.04685735 12.9896112 3.742783e-14 1.646825e-12 A10
## 20 31 10.14299516 13.5602211 3.490344e-13 5.119172e-12 A10
## 21 32 10.04460532 13.3093630 1.626586e-13 3.578490e-12 A10
## 22 33 9.97305515 14.9878564 1.469255e-06 6.464720e-06 A10print(b100.biomass.results)## dap conf.int.low conf.int.up pvalue qvalue genotype
## 1 11 -0.10273614 0.01603354 1.433897e-01 2.175569e-01 B100
## 2 12 -0.05331944 0.14763773 3.417760e-01 3.957407e-01 B100
## 3 13 -0.19433167 0.11252372 5.892521e-01 6.286432e-01 B100
## 4 14 -0.17223071 0.31125464 5.550045e-01 6.105049e-01 B100
## 5 15 -0.56141612 0.03872965 8.427560e-02 1.426203e-01 B100
## 6 16 -0.16839758 1.01001512 1.508052e-01 2.199402e-01 B100
## 7 17 -1.20761685 0.83516601 7.032150e-01 7.195688e-01 B100
## 8 18 -0.16852448 1.69244380 1.035239e-01 1.687056e-01 B100
## 9 19 -1.11131797 1.07162908 9.703722e-01 9.703722e-01 B100
## 10 20 0.95683425 3.18300525 8.908319e-04 2.062979e-03 B100
## 11 21 -1.02609382 3.09523737 3.076201e-01 3.759801e-01 B100
## 12 22 -0.71744932 4.15312089 1.549578e-01 2.199402e-01 B100
## 13 23 -1.72716481 7.54585618 1.698545e-01 2.335500e-01 B100
## 14 25 0.30686251 6.13591628 3.222730e-02 5.672006e-02 B100
## 15 26 -2.16265890 6.18863717 3.194544e-01 3.798917e-01 B100
## 16 27 4.43237220 9.24389988 8.819553e-06 3.185671e-05 B100
## 17 28 1.33865042 8.29025819 9.646199e-03 1.768470e-02 B100
## 18 29 2.22560476 8.60958217 2.396530e-03 4.793060e-03 B100
## 19 30 2.93818406 10.03007759 1.320144e-03 2.904316e-03 B100
## 20 31 8.43451895 12.61765594 5.634832e-10 3.541894e-09 B100
## 21 32 4.53402087 12.71198560 4.007013e-04 1.037109e-03 B100
## 22 33 8.36684530 14.24962413 1.865143e-06 7.460572e-06 B100Output the A10 and B100 tables to the same file
write.table(a10.biomass.results, file='biomass.stats.csv', sep = ',',
row.names = FALSE, append = FALSE)
write.table(b100.biomass.results, file='biomass.stats.csv', sep = ',',
row.names = FALSE, append = TRUE, col.names = FALSE)Biomass response to drought
Calculate the difference in biomass per genotype per day between the 33% and 100% water treatment groups.
drought_response = function(vis.data, genotypes) {
# Initialize drought response data frame with days
drought_resp = data.frame(day=c(
min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))))
# Initialize genotypes
for(g in genotypes) {
control = paste(g,'control',sep='')
drought = paste(g,'drought',sep='')
# Calculate median biomass per treatment per day
drought_resp[,paste(control)] = 0
drought_resp[,paste(drought)] = 0
for(d in min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))) {
drought_resp[drought_resp$day == d, paste(control)] =
median(vis.data[vis.data$genotype == g &
vis.data$treatment == 100 &
as.integer(vis.data$dap) == d,]$fw_biomass)
drought_resp[drought_resp$day == d, paste(drought)] =
median(vis.data[vis.data$genotype == g &
vis.data$treatment == 33 &
as.integer(vis.data$dap) == d,]$fw_biomass)
}
}
drought_resp = drought_resp[-14,]
# Calculate biomass loss to drought
days = c()
response = c()
genotype=c()
for(r in 1:nrow(drought_resp)) {
days = c(days, rep(drought_resp[r,]$day,length(genotypes)))
for(g in genotypes) {
control = paste(g,'control',sep='')
drought = paste(g,'drought',sep='')
response = c(response, drought_resp[r,paste(drought)] -
drought_resp[r,paste(control)])
genotype = c(genotype, g)
}
}
dr.df = data.frame(day=days,response=response,genotype=as.factor(genotype))
return(dr.df)
}A10 vs B100 drought responses
drought.response.parents = drought_response(vis.data.zoom,c('A10','B100'))Plot drought response results
genotype.colors = c('#E6A024','#4AB859')
resp.plot = ggplot(drought.response.parents,
aes(x=day, y=response, group=genotype, color=genotype)) +
geom_point(size=2.5) +
geom_smooth(method="loess") +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(name="Reduced accumulated biomass (g)") +
theme_bw() +
theme(legend.position=c(0.8,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
scale_colour_manual(values=genotype.colors) +
labs(color='Genotype')resp.plot = resp.plot + labs(title='Figure 4C')
print(resp.plot)
Growth curve analysis
In this experiment plants appear to follow relatively basic asymptotic growth patterns. Here we use non-linear least squares regression analysis to model a three-component logistic growth function for each genotype x treatment group.
First load three functions for non-linear growth curve analysis reproduced from:
Paine CET, Marthews TR, Vogt DR, Purves D, Rees M, Hector A, Turnbull LA (2012) How to fit nonlinear plant growth models and calculate growth rates: an update for ecologists. Methods in Ecology and Evolution 3: 245–256. doi: 10.1111/j.2041-210X.2011.00155.x.
output.logis.nlsList <- function(fit, times, CI = F, LOG = F, alpha = 0.05){
coef <- coef(fit)
params <- transform_param.logis(coef)
rates <- list()
groups <- rownames(params)
n.groups <- nrow(coef)
# compute rates for each group seperately
rates <- list()
for(i in 1:(n.groups)){
K <- params[i,1]; r <- params[i,2]; M0 <- params[i,3]
rates[[i]] = data.frame(
times = times,
M = (M0*K)/(M0+(K-M0)*exp(-r*times)),
AGR = (r*M0*K*(K-M0)*exp(-r*times))/(M0+(K-M0)*exp(-r*times))^2
)
rates[[i]]$RGRt <- rates[[i]]$AGR/rates[[i]]$M
rates[[i]]$RGRm <- r*(1 - rates[[i]]$M/K)
if(LOG == T){
rates[[i]]$RGRt <- rates[[i]]$AGR
rates[[i]]$RGRm <- r*rates[[i]]$M*(1-rates[[i]]$M/K)
rates[[i]]$AGR <- rates[[i]]$AGR*exp(rates[[i]]$M)
}
# commute CIs for each group's estaimates, if desired
if(CI == T){
cov <- summary(fit)$cov[i,,]
x <- y <- data.frame(rmvnorm(n=1000, mean=c(coef[i, 1], coef[i, 2], coef[i, 3]), sigma=cov))
x$K <- y[,1]
x$r <- 1/y[,3]
x$M0 <- y[,1]/(1 + exp(y[,2]/y[,3]))
M <- AGR <- RGRt <- RGRm <- matrix(NA, ncol = length(times), nrow = nrow(x))
for(j in 1:nrow(x)){
K <- x[j,4]; r <- x[j,5]; M0 <- x[j,6]
M[j,] <- (M0*K)/(M0+(K-M0)*exp(-r*times))
AGR[j,] <- (r*M0*K*(K-M0)*exp(-r*times))/(M0+(K-M0)*exp(-r*times))^2
RGRt[j,] <- AGR[j,]/M[j,]
RGRm[j,] <- r*(1 - M[j,]/K)
if(LOG ==T){
RGRt[j,] <- AGR[j,]
RGRm[j,] <- r*M[j,]*(1 - M[j,]/K)
AGR[j,] <- AGR[j,]*exp(M[j,])
}
}
}
CIs <- summarizer(list(M = M, AGR = AGR, RGRt = RGRt, RGRm = RGRm), alpha)
rates[[i]] <- cbind(rates[[i]], CIs)
}
names(rates) <- rownames(params)
# now compute differences among groups
diffs <- list()
for(i in 1:(n.groups-1)){
Ki <- params[i,1]; ri <- params[i,2]; M0i <- params[i,3]
for(j in (i+1):n.groups){
Kj <- params[j,1]; rj <- params[j,2]; M0j <- params[j,3]
diffs.ij = data.frame(
times = times,
diffM = rates[[i]]$M - rates[[j]]$M,
diffAGR = rates[[i]]$AGR - rates[[j]]$AGR,
diffRGRt = rates[[i]]$RGRt - rates[[j]]$RGRt
)
# comparing RGRm has to be done on a biomass basis. So it needs special treatment. First, we aev to know what range of biomasses are shared between groups.
Mmin <- max(min(rates[[i]]$M), min(rates[[j]]$M)) # yieds the range of overlapping masses between groups i & j
Mmax <- min(max(rates[[i]]$M), max(rates[[j]]$M))
#diffs.ij$Mseq <- seq(Mmin, Mmax, length = 100)
diffs.ij$Mseq <- seq(Mmin, Mmax, length = length(times))
if(LOG == F){
diffs.ij$diffRGRm <- ri*(1 - diffs.ij$Mseq/Ki) - rj*(1 - diffs.ij$Mseq/Kj)
} else{
diffs.ij$diffRGRm <- ri*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Ki) - rj*Mseq*(1 - diffs.ij$Mseq/Kj)
}
}
if(CI == T){
# get params for group i
covi <- summary(fit)$cov[i,,]
xi <- yi <- data.frame(rmvnorm(n=1000, mean=c(coef[i, 1], coef[i, 2], coef[i, 3]), sigma=covi))
xi$K <- yi[,1]
xi$r <- 1/yi[,3]
xi$M0 <- yi[,1]/(1 + exp(yi[,2]/yi[,3]))
# get params for group j
covj <- summary(fit)$cov[j,,]
xj <- yj <- data.frame(rmvnorm(n=1000, mean=c(coef[j, 1], coef[j, 2], coef[j, 3]), sigma=covj))
xj$K <- yj[,1]
xj$r <- 1/yj[,3]
xj$M0 <- yj[,1]/(1 + exp(yj[,2]/yj[,3]))
# now compute diffs for each random set of drawn parameters
Mi <- Mj <- AGRi <- AGRj <- RGRti <- RGRtj <- RGRmi <- RGRmj <- diffM <- diffAGR <- diffRGRt <- diffRGRm <- matrix(NA, ncol = length(times), nrow = nrow(xi))
for(k in 1:nrow(xi)){
Ki <- xi[k,4]; ri <- xi[k,5]; M0i <- xi[k,6]
Kj <- xj[k,4]; rj <- xj[k,5]; M0j <- xj[k,6]
Mi[k,] <- (M0i*Ki)/(M0i+(Ki-M0i)*exp(-ri*times))
Mj[k,] <- (M0j*Kj)/(M0j+(Kj-M0j)*exp(-rj*times))
AGRi[k,] <- (ri*M0i*Ki*(Ki-M0i)*exp(-ri*times))/(M0i+(Ki-M0i)*exp(-ri*times))^2
AGRj[k,] <- (rj*M0j*Kj*(Kj-M0j)*exp(-rj*times))/(M0j+(Kj-M0j)*exp(-rj*times))^2
RGRti[k,] <- AGRi[k,]/Mi[k,]
RGRtj[k,] <- AGRj[k,]/Mj[k,]
RGRmi[k,] <- ri*(1 - diffs.ij$Mseq/Ki)
RGRmj[k,] <- rj*(1 - diffs.ij$Mseq/Kj)
if(LOG == T){
RGRti[k,] <- AGRi[k,]
RGRtj[k,] <- AGRj[k,]
RGRmi[k,] <- ri*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Ki)
RGRmj[k,] <- rj*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Kj)
AGRi[k,] <- AGRi[k,]*exp(Mi[k,])
AGRj[k,] <- AGRj[k,]*exp(Mj[k,])
}
diffM[k,] <- Mi[k,] - Mj[k,]
diffAGR[k,] <- AGRi[k,] - AGRj[k,]
diffRGRt[k,] <- RGRti[k,] - RGRtj[k,]
diffRGRm[k,] <- RGRmi[k,] - RGRmj[k,]
}
CIs <- summarizer(list(diffM = diffM, diffAGR = diffAGR, diffRGRt = diffRGRt, diffRGRm = diffRGRm), alpha)
diffs[[paste(groups[i], groups[j], sep = "_")]] <- cbind(diffs.ij, CIs)
} else{
diffs[[paste(groups[i], groups[j], sep = "_")]] <- diffs.ij
}
} # end loop over pairwise combinations of groups
out <- list(params = params, rates = rates, diffs = diffs)
return(out)
}
transform_param.logis <- function(coef){
K = coef[1]
r = 1/(coef[3])
M0 = K/(1 + exp(coef[2]/coef[3])) #untransform best-fit parameters to K, r and M0
if(is.data.frame(K)){
out <- cbind(K, r, M0)
} else {
out <- c(K, r, M0)
}
names(out) <- c("K", "r", "M0")
return(out)
}
# this function returns confidence envelopes around growth trajectories, and growth rates.
summarizer <- function(dat, alpha){
n <- length(dat)
quantiles <- c(alpha/2, 1-(alpha/2))
CIs <- data.frame(matrix(NA, ncol(dat[[1]]), n*2))
names(CIs) <- paste(rep(names(dat), each = 2), c("lo", "hi"), sep = ".")
for(i in 1:n){
CIs[,(2*i-1):(2*i)] <- t(apply(dat[[i]], 2, quantile, quantiles, na.rm = T))
}
return(CIs)
}Subset the VIS data
sub.vis.data = vis.data.zoom[(vis.data.zoom$treatment == 100 | vis.data.zoom$treatment == 33),
c("plant_id","dap","fw_biomass","group")]
sub.vis.data$group = factor(sub.vis.data$group)Group data
grouped.sub.vis.data = groupedData(fw_biomass ~ dap | group, sub.vis.data)Fit three-component logistic functions for each genotype x treatment group
fit.logis = nlsList(fw_biomass ~ SSlogis(dap, Asym, xmid, scal),
data = grouped.sub.vis.data)Output estimated growth rate parameters at even time intervals
est.interval = seq(min(sub.vis.data$dap), max(sub.vis.data$dap), length = 100)
out.fit.logis = output.logis.nlsList(fit.logis, times = est.interval,
CI = TRUE, LOG = FALSE, alpha = 0.05)Compute the time and magnitude of maximum growth rate
logis.results = data.frame(group=levels(sub.vis.data$group))
logis.results$max.AGR = 0
logis.results$max.AGR.dap = 0
logis.results$AGR.lo = 0
logis.results$AGR.hi = 0
group.order = names(out.fit.logis$rates)
for(i in 1:length(group.order)) {
logis.results[logis.results == group.order[i],]$max.AGR.dap =
out.fit.logis$rates[[i]]$times[out.fit.logis$rates[[i]]$AGR ==
max(out.fit.logis$rates[[i]]$AGR)]
logis.results[logis.results == group.order[i],]$max.AGR =
max(out.fit.logis$rates[[i]]$AGR)
logis.results[logis.results == group.order[i],]$AGR.lo =
out.fit.logis$rates[[i]]$AGR.lo[out.fit.logis$rates[[i]]$AGR ==
max(out.fit.logis$rates[[i]]$AGR)]
logis.results[logis.results == group.order[i],]$AGR.hi =
out.fit.logis$rates[[i]]$AGR.hi[out.fit.logis$rates[[i]]$AGR ==
max(out.fit.logis$rates[[i]]$AGR)]
}
print(logis.results)## group max.AGR max.AGR.dap AGR.lo AGR.hi
## 1 A10-100 3.031821 24.23556 2.962031 3.099600
## 2 A10-33 1.800921 25.36054 1.724413 1.889975
## 3 B100-100 2.428901 26.48553 2.378344 2.479694
## 4 B100-33 1.693055 24.68555 1.604441 1.784634
## 5 R102-100 2.679992 28.06050 2.596513 2.757770
## 6 R102-33 1.890607 26.48553 1.813447 1.977527
## 7 R128-100 3.090714 25.58554 3.003102 3.174770
## 8 R128-33 1.847410 24.46055 1.770401 1.930073
## 9 R133-100 3.380220 23.11057 3.222844 3.526918
## 10 R133-33 2.292966 23.33557 2.180286 2.421985
## 11 R161-100 2.952080 25.58554 2.853000 3.043931
## 12 R161-33 1.813290 24.68555 1.729947 1.913824
## 13 R187-100 2.535501 24.46055 2.451728 2.629045
## 14 R187-33 1.676044 24.46055 1.596170 1.763532
## 15 R20-100 3.360632 22.88557 3.225079 3.499300
## 16 R20-33 1.949526 22.66058 1.858486 2.053219
## 17 R70-100 2.978971 24.01056 2.886229 3.083468
## 18 R70-33 1.654110 23.56057 1.571712 1.742591
## 19 R98-100 2.597831 24.23556 2.507745 2.691376
## 20 R98-33 1.538670 24.68555 1.467270 1.612429Plot logistic growth models for S. viridis and S. italica and absolute growth rates over time
parent.groups = c("A10-100", "A10-33", "B100-100", "B100-33")
group.colors = c("#F8766D", "#7CAE00", "#00BFC4", "#C77CFF")
# Biomass over time
gca.plot = ggplot(vis.data.zoom[(vis.data.zoom$genotype == 'A10' |
vis.data.zoom$genotype == 'B100') &
(vis.data.zoom$treatment == 100 | vis.data.zoom$treatment == 33),],
aes(x=dap, y=fw_biomass, color=factor(group))) +
geom_point(size=2.5) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-1,41), name="Estimated biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
scale_colour_manual(name='Genotype-treatment',
values=group.colors, labels=parent.groups)
agr.plot = ggplot() +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(name="Absolute growth rate (g/day)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
scale_colour_manual(name='Genotype-treatment',
values=group.colors, labels=parent.groups)
# Add logistic growth models and confidence intervals
for(g in 1:length(parent.groups)) {
group = parent.groups[g]
i = match(group, group.order)
group.rates = as.data.frame(out.fit.logis$rates[i])
col.name = gsub('-', '.', group)
conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
rev(group.rates[,paste(col.name, '.times', sep='')])),
y=c(group.rates[,paste(col.name, '.M.lo', sep='')],
rev(group.rates[,paste(col.name, '.M.hi', sep='')])))
agr.conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
rev(group.rates[,paste(col.name, '.times', sep='')])),
y=c(group.rates[,paste(col.name, '.AGR.lo', sep='')],
rev(group.rates[,paste(col.name, '.AGR.hi', sep='')])))
gca.plot = gca.plot + geom_polygon(data=conf.int, aes(x=x, y=y),
fill='gray60', color=NA, alpha=0.4) +
geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
y=paste(col.name, '.M', sep='')),
color=group.colors[g])
agr.plot = agr.plot +
geom_polygon(data=agr.conf.int, aes(x=x, y=y),
fill='gray60', color=NA, alpha=0.4) +
geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
y=paste(col.name, '.AGR', sep='')),
color=group.colors[g])
}gca.plot = gca.plot + labs(title='Figure 4B')
print(gca.plot)
agr.plot = agr.plot + labs(title='Figure 4C')
print(agr.plot)
Plot logistic growth models for all Setaria genotypes
plot_gca = function(genotype, vis.data, out.fit.logis) {
groups = c(paste(genotype,'-100',sep=''), paste(genotype,'-33',sep=''))
group.colors = c("#00BFC4", "#F8766D")
# Biomass over time
gca.plot = ggplot(vis.data[(vis.data$genotype == genotype) &
(vis.data$treatment == 100 | vis.data$treatment == 33),],
aes(x=dap, y=fw_biomass, color=factor(treatment))) +
geom_point(size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-1,41), name="Estimated biomass (g)") +
theme_bw() +
theme(legend.position='none',
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Treatment") +
facet_grid(. ~ genotype)
# Add logistic growth models and confidence intervals
for(g in 1:length(groups)) {
group = groups[g]
i = match(group, group.order)
group.rates = as.data.frame(out.fit.logis$rates[i])
col.name = gsub('-', '.', group)
conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
rev(group.rates[,paste(col.name, '.times', sep='')])),
y=c(group.rates[,paste(col.name, '.M.lo', sep='')],
rev(group.rates[,paste(col.name, '.M.hi', sep='')])))
gca.plot = gca.plot + geom_polygon(data=conf.int, aes(x=x, y=y),
fill='gray60', color=NA, alpha=0.4) +
geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
y=paste(col.name, '.M', sep='')),
color=group.colors[g])
}
return(gca.plot)
}Multiple plot function modified from the Cookbook for R by Winston Chang: http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/
multiplot <- function(..., plotlist=NULL, cols=1, layout=NULL) {
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols),byrow=TRUE)
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}Generate the Setaria biomass/growth plots
a10.gca.plot = plot_gca('A10', vis.data.zoom, out.fit.logis)
b100.gca.plot = plot_gca('B100', vis.data.zoom, out.fit.logis)
r102.gca.plot = plot_gca('R102', vis.data.zoom, out.fit.logis)
r128.gca.plot = plot_gca('R128', vis.data.zoom, out.fit.logis)
r133.gca.plot = plot_gca('R133', vis.data.zoom, out.fit.logis)
r161.gca.plot = plot_gca('R161', vis.data.zoom, out.fit.logis)
r187.gca.plot = plot_gca('R187', vis.data.zoom, out.fit.logis)
r20.gca.plot = plot_gca('R20', vis.data.zoom, out.fit.logis)
r70.gca.plot = plot_gca('R70', vis.data.zoom, out.fit.logis)
r98.gca.plot = plot_gca('R98', vis.data.zoom, out.fit.logis)Supplemental Figure S3
multiplot(a10.gca.plot, b100.gca.plot, r102.gca.plot, r128.gca.plot,
r133.gca.plot, r161.gca.plot, r187.gca.plot, r20.gca.plot,
r70.gca.plot, r98.gca.plot, cols=3)
agr.plot = agr.plot + labs(title='Figure 4C')
print(agr.plot)
Water use efficiency
In this section we combine data on the volume of water added to each plant per water application (up to the target weight) with fresh-weight biomass data to estimate water use efficiency.
Download water volume data (n = 33,496 water applications).
if (!file.exists('B2_watering_data_phenofront.csv')) {
download.file('http://files.figshare.com/1845363/B2_watering_data_phenofront.csv',
'B2_watering_data_phenofront.csv')
}Read data and format for analysis
water.data = read.table(file="B2_watering_data_phenofront.csv", sep=",", header=TRUE)Remove the empty pot controls (Barcodes start with zero).
water.data = water.data[grep('000A', water.data$plant.barcode, invert = TRUE),]Remove snapshots that were not valid waterings (water.amount = -1).
water.data = water.data[water.data$water.amount != -1,]Add a treatment group column coded in the barcode.
water.data$treatment <- NA
water.data$treatment[grep("AA", water.data$plant.barcode)] <- 100
water.data$treatment[grep("AB", water.data$plant.barcode)] <- 0
water.data$treatment[grep("AC", water.data$plant.barcode)] <- 16
water.data$treatment[grep("AD", water.data$plant.barcode)] <- 33
water.data$treatment[grep("AE", water.data$plant.barcode)] <- 66Add a plant genotype column coded in the barcode.
water.data$genotype <- NA
water.data$genotype[grep("p1", water.data$plant.barcode)] <- 'A10'
water.data$genotype[grep("p2", water.data$plant.barcode)] <- 'B100'
water.data$genotype[grep("r1", water.data$plant.barcode)] <- 'R20'
water.data$genotype[grep("r2", water.data$plant.barcode)] <- 'R70'
water.data$genotype[grep("r3", water.data$plant.barcode)] <- 'R98'
water.data$genotype[grep("r4", water.data$plant.barcode)] <- 'R102'
water.data$genotype[grep("r5", water.data$plant.barcode)] <- 'R128'
water.data$genotype[grep("r6", water.data$plant.barcode)] <- 'R133'
water.data$genotype[grep("r7", water.data$plant.barcode)] <- 'R161'
water.data$genotype[grep("r8", water.data$plant.barcode)] <- 'R187'Add a genotype x treatment group column.
water.data$group = paste(water.data$genotype,'-',water.data$treatment,sep='')Encode the calendar time column as a date-time.
water.data$timestamp = ymd_hms(water.data$timestamp)Add a column for days after planting since the planting date.
water.data$dap = as.numeric(water.data$timestamp - planting_date)Some plants were removed from the system but were not inactivated, so remove them.
bad.cars = unique(water.data$car.tag[water.data$water.amount >120 &
water.data$dap > 14])
water.data = water.data[!water.data$car.tag %in% bad.cars,]Data for Figure 1B
Extract water data for S. viridis for Figure 1B.
sv.water = water.data[(water.data$treatment == 100 | water.data$treatment == 33) &
water.data$genotype == 'A10',]Classify water job as early or later in the day. Early water treatments started at ZT23, later treatments started at ZT8.
sv.water$early.late <- NA #restrict to the scheduled watering later in the run
sv.water$early.late[hour(sv.water$timestamp) < 7] <- 'ZT23'
sv.water$early.late[hour(sv.water$timestamp) < 15 &
hour(sv.water$timestamp) > 11 ] <- 'ZT8'
sv.water = sv.water[!is.na(sv.water$early.late),]
sv.water$dap = day(sv.water$timestamp) + 4
sv.water = sv.water[sv.water$dap > 14,]Plot the S. viridis water volume data
water.plot = ggplot(sv.water, aes(y = water.amount, x = dap,
group = factor(interaction(early.late,
treatment,dap)),
fill = factor(interaction(early.late,
treatment)),
color = factor(interaction(early.late,
treatment)))) +
geom_boxplot(outlier.colour = NULL) +
scale_x_continuous("Days after planting") +
scale_y_continuous("Water volume (ml)") +
theme_bw() +
theme(legend.position = c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(fill='Water treatment', color='Water treatment')water.plot = water.plot + labs(title='Figure 1B')
print(water.plot)
WUE function
Calculate WUE for S. viridis and S. italica
wue.parents = function(water, biomass) {
# WUE data vectors
dap.list = c()
plant.list = c()
water.list = c()
biomass.list = c()
treatment.list = c()
group.list=c()
# Get unique barcodes for biomass
barcodes = unique(biomass[(biomass$genotype=='A10' | biomass$genotype=='B100') &
(biomass$treatment == 100 |
biomass$treatment == 33),]$plant_id)
for(barcode in barcodes) {
snapshots = biomass[biomass$plant_id==barcode,]
snapshots = snapshots[with(snapshots, order(dap)),]
for(row in 1:nrow(snapshots)) {
total.water = sum(water[water$dap <= snapshots[row,]$dap &
water$plant.barcode == barcode,]$water.amount)
dap.list = c(dap.list, snapshots[row,]$dap)
plant.list = c(plant.list, barcode)
water.list = c(water.list, total.water)
biomass.list = c(biomass.list, snapshots[row,]$fw_biomass)
treatment.list = c(treatment.list, snapshots[row,]$treatment)
group.list = c(group.list,snapshots[row,]$group)
}
}
wue.data = data.frame(plant_id=plant.list,
dap=dap.list,
water=water.list,
biomass=biomass.list,
treatment=treatment.list,
group=group.list)
return(wue.data)
}
parents.wue = wue.parents(water.data, vis.data.zoom)
parents.wue = parents.wue[parents.wue$water != 0,]Plot S. viridis and S. italica WUE in mg/mL.
wue.plot = ggplot(parents.wue, aes(x=dap, y=(biomass/water)*1000,
color=factor(group))) +
geom_point(size=2.5) +
geom_smooth(method="loess", size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-2.1, 26),
name="Water-use efficiency (mg/mL)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color='Genotype-treatment')wue.plot = wue.plot + labs(title='Figure 4D')
print(wue.plot)
Statistical analysis of WUE
Genotype WUE function
wue = function(water, biomass, genotype) {
# WUE data vectors
dap.list = c()
plant.list = c()
water.list = c()
biomass.list = c()
treatment.list = c()
# Get unique barcodes for biomass
barcodes = unique(biomass[biomass$genotype==genotype &
(biomass$treatment == 100 |
biomass$treatment == 33),]$plant_id)
for(barcode in barcodes) {
snapshots = biomass[biomass$plant_id==barcode,]
snapshots = snapshots[with(snapshots, order(dap)),]
for(row in 1:nrow(snapshots)) {
total.water = sum(water[water$dap <= snapshots[row,]$dap &
water$plant.barcode == barcode,]$water.amount)
dap.list = c(dap.list, snapshots[row,]$dap)
plant.list = c(plant.list, barcode)
water.list = c(water.list, total.water)
biomass.list = c(biomass.list, snapshots[row,]$fw_biomass)
treatment.list = c(treatment.list, snapshots[row,]$treatment)
}
}
wue.data = data.frame(plant_id=plant.list,
dap=dap.list,
water=water.list,
biomass=biomass.list,
treatment=treatment.list)
return(wue.data)
}Analyse pair-wise treatment differences for two-day increments for each genotype.
analyze_wue = function(wue.data) {
days = c()
diff.low = c()
diff.up = c()
pvals = c()
for(day in levels(factor(as.integer(wue.data$dap)))) {
day = as.integer(day)
control = wue.data[(as.integer(wue.data$dap) == day |
as.integer(wue.data$dap) == day + 1) &
wue.data$treatment == 100,]
drought = wue.data[(as.integer(wue.data$dap) == day |
as.integer(wue.data$dap) == day + 1) &
wue.data$treatment == 33,]
control.wue = control$biomass / control$water
drought.wue = drought$biomass / drought$water
test = t.test(x=control.wue,y=drought.wue)
days = c(days,day)
diff.low = c(diff.low, test$conf.int[1])
diff.up = c(diff.up, test$conf.int[2])
pvals = c(pvals, test$p.value)
}
results = data.frame(dap=as.numeric(days),
conf.int.low=diff.low,
conf.int.up=diff.up,
pvalue=pvals)
return(results)
}Analyze WUE for each genotype and aggregate results
wue.results = data.frame()
for(genotype in c('A10', 'B100')) {
genotype.wue = wue(water.data, vis.data.zoom, genotype)
genotype.wue = genotype.wue[genotype.wue$water != 0,]
genotype.wue.results = analyze_wue(genotype.wue)
genotype.wue.results$genotype = genotype
wue.results = rbind(wue.results, genotype.wue.results)
}Control for multiple testing by controlling the FDR
qvalues.wue = p.adjust(wue.results$pvalue,method="fdr")
wue.results$qvalue = qvalues.wue
write.table(wue.results, file='wue.stats.csv', sep=',', row.names=FALSE)
print(wue.results)## dap conf.int.low conf.int.up pvalue genotype qvalue
## 1 11 0.0001681857 7.003870e-04 0.0023589893 A10 0.03459851
## 2 12 -0.0002932779 6.855075e-04 0.4161121704 A10 0.79604067
## 3 13 -0.0003981046 7.310925e-04 0.5504252348 A10 0.85019455
## 4 14 -0.0003666556 1.086597e-03 0.3151858236 A10 0.65930472
## 5 15 -0.0006318583 1.290448e-03 0.4824013893 A10 0.85019455
## 6 16 -0.0011579704 1.505706e-03 0.7841660738 A10 0.85019455
## 7 17 -0.0030117768 2.033894e-03 0.6757645778 A10 0.85019455
## 8 18 -0.0019852804 1.537689e-03 0.7920861585 A10 0.85019455
## 9 19 -0.0022189870 1.423300e-03 0.6543284916 A10 0.85019455
## 10 20 -0.0018767476 2.487869e-03 0.7709120024 A10 0.85019455
## 11 21 -0.0011815956 3.333287e-03 0.3296523601 A10 0.65930472
## 12 22 -0.0011941099 3.447888e-03 0.3180544763 A10 0.65930472
## 13 23 -0.0030785137 4.736555e-03 0.6313454757 A10 0.85019455
## 14 25 -0.0017208810 2.529517e-03 0.6913194980 A10 0.85019455
## 15 26 -0.0016761719 2.370431e-03 0.7214693240 A10 0.85019455
## 16 27 -0.0015099451 1.887345e-03 0.8191576027 A10 0.85816511
## 17 28 -0.0016265463 1.712963e-03 0.9574064089 A10 0.95740641
## 18 29 -0.0016959640 1.418049e-03 0.8545599052 A10 0.87443339
## 19 30 -0.0017780029 1.240136e-03 0.7142953563 A10 0.85019455
## 20 31 -0.0018270439 1.050066e-03 0.5808623905 A10 0.85019455
## 21 32 -0.0018669052 9.419976e-04 0.5014284454 A10 0.85019455
## 22 33 -0.0019935110 2.547949e-03 0.7922267419 A10 0.85019455
## 23 11 -0.0004686284 1.653400e-04 0.3231119733 B100 0.65930472
## 24 12 -0.0003176706 4.659287e-04 0.6970305202 B100 0.85019455
## 25 13 -0.0006793007 2.052254e-04 0.2809421281 B100 0.65930472
## 26 14 -0.0005693316 7.964500e-04 0.7298612911 B100 0.85019455
## 27 15 -0.0013429980 2.693045e-04 0.1789099092 B100 0.56228829
## 28 16 -0.0011018089 1.726678e-03 0.6393283914 B100 0.85019455
## 29 17 -0.0053976734 -4.565038e-04 0.0249473645 B100 0.11089327
## 30 18 -0.0036398672 2.934236e-04 0.0894445265 B100 0.30273532
## 31 19 -0.0050384129 -6.422642e-04 0.0147704390 B100 0.08123741
## 32 20 -0.0044244661 -5.891914e-05 0.0447774175 B100 0.17910967
## 33 21 -0.0064021316 -1.321929e-03 0.0053896550 B100 0.03952414
## 34 22 -0.0057466915 -1.176396e-03 0.0051207769 B100 0.03952414
## 35 23 -0.0062955680 -1.663831e-03 0.0033048057 B100 0.03635286
## 36 25 -0.0058026119 -1.706030e-03 0.0009451266 B100 0.03381652
## 37 26 -0.0052016475 -1.378281e-03 0.0015371147 B100 0.03381652
## 38 27 -0.0043036325 -5.750897e-04 0.0125215176 B100 0.07870668
## 39 28 -0.0037893424 -2.728566e-04 0.0252030152 B100 0.11089327
## 40 29 -0.0033915507 1.520038e-04 0.0711932892 B100 0.26104206
## 41 30 -0.0027631015 6.920535e-04 0.2276453890 B100 0.65930472
## 42 31 -0.0025651583 7.463200e-04 0.2647434581 B100 0.65930472
## 43 32 -0.0025173692 8.428216e-04 0.3102535812 B100 0.65930472
## 44 33 -0.0029600383 2.256215e-03 0.7752593144 B100 0.85019455Model tiller number
In this section we use linear modeling to estimate tiller number.
Calculate height-width ratio.
vis.data.zoom$height_width_ratio = vis.data.zoom$height_above_bound / vis.data.zoom$extent_xDownload data for manually measured tiller counts 195 random images.
if (!file.exists('200tiller_counts.txt')) {
download.file('http://files.figshare.com/1845359/200tiller_counts.txt',
'200tiller_counts.txt')
}Read manual tiller count data.
tiller.counts = read.table(file="200tiller_counts.txt",sep="\t",header=TRUE)Format date.
tiller.counts$date = as.POSIXct(paste(paste(tiller.counts$year,
tiller.counts$month,
tiller.counts$day, sep='-'),
paste(tiller.counts$hours,
tiller.counts$minutes,
tiller.counts$seconds, sep=':'), sep=' '))Merge with VIS snapshot table.
tiller.counts = merge(x=tiller.counts, y=vis.data.zoom, by = 'date')Model tiller counts.
tiller.model.full = lm(ave_tillers ~ fw_biomass + height_above_bound + solidity +
extent_x + height_width_ratio, data=tiller.counts)Variance inflation factor.
vif(tiller.model.full)## fw_biomass height_above_bound solidity
## 16.615475 21.832350 1.591739
## extent_x height_width_ratio
## 14.738418 3.238483Drop height and recalculate VIF
tiller.model1 = lm(ave_tillers ~ fw_biomass + solidity + extent_x +
height_width_ratio, data=tiller.counts)
vif(tiller.model1)## fw_biomass solidity extent_x
## 10.014163 1.363769 11.118321
## height_width_ratio
## 1.847957Drop extent x.
tiller.model2 = lm(ave_tillers ~ fw_biomass + solidity + height_width_ratio,
data=tiller.counts)
vif(tiller.model2)## fw_biomass solidity height_width_ratio
## 1.024829 1.041759 1.039990summary(tiller.model2)##
## Call:
## lm(formula = ave_tillers ~ fw_biomass + solidity + height_width_ratio,
## data = tiller.counts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6877 -1.2636 -0.0169 0.8454 5.6053
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.23276 0.54698 9.567 < 2e-16 ***
## fw_biomass 0.21998 0.01254 17.538 < 2e-16 ***
## solidity 0.15006 1.91657 0.078 0.938
## height_width_ratio -2.19370 0.32640 -6.721 2.12e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.707 on 187 degrees of freedom
## Multiple R-squared: 0.644, Adjusted R-squared: 0.6382
## F-statistic: 112.7 on 3 and 187 DF, p-value: < 2.2e-16Drop solidity
tiller.model3 = lm(ave_tillers ~ fw_biomass + height_width_ratio, data=tiller.counts)
summary(tiller.model3)##
## Call:
## lm(formula = ave_tillers ~ fw_biomass + height_width_ratio, data = tiller.counts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.6886 -1.2680 -0.0274 0.8433 5.6022
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.26031 0.41771 12.593 < 2e-16 ***
## fw_biomass 0.21986 0.01242 17.709 < 2e-16 ***
## height_width_ratio -2.18932 0.32072 -6.826 1.17e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.703 on 188 degrees of freedom
## Multiple R-squared: 0.6439, Adjusted R-squared: 0.6402
## F-statistic: 170 on 2 and 188 DF, p-value: < 2.2e-16Plot partial regressions.
avPlots(model = tiller.model3, pch=16, layout=c(2,1), main="Figure 5B")
Download data for manually measured tiller counts for 646 random images.
if (!file.exists('660tiller_counts.txt')) {
download.file('http://files.figshare.com/1845358/660tiller_counts.txt',
'660tiller_counts.txt')
}Predict tiller count for a second set of ~660 randomly selected plants.
tiller660 = read.table(file="660tiller_counts.txt", sep="\t", header=TRUE)
tiller660 = merge(x=tiller660, y=vis.data.zoom, by='datetime')
tiller660$predicted_tiller_count = predict.lm(object = tiller.model3,
newdata=tiller660)Summarize the difference between the predicted and manual tiller counts
tiller.diff = tiller660$predicted_tiller_count - tiller660$tiller_count
t.test(tiller.diff)##
## One Sample t-test
##
## data: tiller.diff
## t = 14.7453, df = 645, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.7978015 1.0429342
## sample estimates:
## mean of x
## 0.9203678Plot manual versus predicted tiller counts.
tiller.plot = ggplot(data=tiller660, aes(x=tiller_count, y=predicted_tiller_count,
color=factor(genotype))) +
geom_point(size=2.5) +
geom_abline(intercept=0, slope=1) +
scale_x_continuous(lim=c(-1,14), "Manual tiller count") +
scale_y_continuous(lim=c(-1,14), "Predicted tiller count") +
theme_bw() +
theme(legend.position=c(0.45,0.95),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"),
legend.direction="horizontal") +
labs(color="Genotype")tiller.plot = tiller.plot + labs(title='Figure 5C')
print(tiller.plot)
Predict tiller counts for the whole data set.
vis.data.zoom$tiller_count = predict.lm(object = tiller.model3, newdata = vis.data.zoom)Plot tiller counts for S. viridis and S. italica water treatments 100% and 33% full-capacity.
tc.par.plot = ggplot(vis.data.zoom[(vis.data.zoom$genotype == 'A10' |
vis.data.zoom$genotype == 'B100') &
(vis.data.zoom$treatment == 100 |
vis.data.zoom$treatment == 33),],
aes(x=dap, y=tiller_count, color=factor(group))) +
geom_point(size=2.5) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(0,12),
name="Estimated tiller count") +
theme_bw() +
theme(legend.position=c(0.25,0.75),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype-treatment")## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 2 rows containing missing values
## (geom_point).tc.par.plot = tc.par.plot + labs(title='Figure 5D')
print(tc.par.plot)## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 2 rows containing missing values
## (geom_point).
Plot estimated tiller counts for all genotypes.
tiller.facet.plot = ggplot(vis.data.zoom[vis.data.zoom$treatment == 100 |
vis.data.zoom$treatment == 33,],
aes(x=dap, y=tiller_count, color=factor(treatment))) +
geom_point(size=1) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(0,12),
name="Estimated tiller count") +
theme_bw() +
theme(legend.position=c(0.8,0.1),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Treatment") +
facet_wrap(~ genotype, ncol = 3)## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 3 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 2 rows containing missing values
## (geom_point).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).## Warning in loop_apply(n, do.ply): Removed 4 rows containing missing values
## (geom_point).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).tiller.facet.plot = tiller.facet.plot + labs(title='Supplemental Figure S4')
print(tiller.facet.plot)## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 3 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (stat_smooth).## Warning in loop_apply(n, do.ply): Removed 2 rows containing missing values
## (geom_point).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).## Warning in loop_apply(n, do.ply): Removed 4 rows containing missing values
## (geom_point).## Warning in loop_apply(n, do.ply): Removed 1 rows containing missing values
## (geom_point).
Trait table
Output a table of all of the traits input and added here. Export all of the VIS traits (except color).
h.table = vis.data.zoom
h.table$plant_id = as.character(h.table$plant_id)
h.table$wue = NA
for(r in 1:nrow(h.table)) {
row = h.table[r,]
total.water = sum(water.data[water.data$dap <= row$dap & water.data$plant.barcode == row$plant_id,]$water.amount)
h.table[r,]$wue = row$sv_area / total.water
}
write.table(h.table,file="vis.traits.csv",quote=FALSE,sep=",",row.names=FALSE)